library(tidysynth)     #Running synthetic control models
library(tidyverse)     #Load tidyverse for piping and ggplot edits

load("Data/synth_data_clean.RData")

#### Synthetic Control Methods at the State-Level ######

# Vaccination rates. DV: Administered_Dose1_Recip_18PlusPop_Pct--Percent of 18+ population with at least one dose based on the jurisdiction where recipient lives, Series_Complete_18PlusPop_Pct -- Percent of people 18+ who are fully vaccinated (have second dose of a two-dose vaccine or one dose of a single-dose vaccine) based on the jurisdiction where recipient lives. Series_Complete_18Plus -- Total number of people 18+ who are fully vaccinated (have second dose of a two-dose vaccine or one dose of a single-dose vaccine) based on the jurisdiction where recipient lives.


###################ANALYSIS###################

### Running 1st dose synthetic control ###

#Creating empty lists to hold results
administered_dose1_trends <- list()
administered_dose1_differences <- list()
administered_dose1_weights <- list()
administered_dose1_balance <- list()
administered_dose1_placebos <- list()
administered_dose1_dat <- list()

#setwd() #Set a working directory for images from the analyses to be saved
for(i in unique(lotteries$state[lotteries$lottery_incentive %in% 1])){
  
  # Take out all untreated states 
  
  x <- subset(state_vaccines,state_vaccines$state %in% i)
  x <- rbind(x,untreated_states)
  x <- x[order(x$date,x$state),]
  
  #take out observations before April 19th (widely available vaccine date)
  x <- x[x$date >= as.Date("2021-04-19"),]
  
  state_vaccines_out <-
    x %>%
    # initial the synthetic control object
    synthetic_control(outcome = Administered_Dose1_Recip_18PlusPop_Pct, # outcome
                      unit = state, # unit index in the panel data
                      time = date, # time index in the panel data
                      i_unit = i, # unit where the intervention occurred
                      i_time = unique(state_vaccines$date[state_vaccines$date %in% as.Date(lotteries$lottery_announced[lotteries$state %in% i])]), # time period when the intervention occurred
                      generate_placebos=T # generate placebo synthetic controls (for inference)
    ) %>%
    
    # Generate the aggregate predictors used to fit the weights
    
    generate_predictor(time_window = as.Date("2021-04-19"),
                       pres_dem_two_party_vote_percent = pres_dem_two_party_vote_percent,
                       income_per_capita=income_per_capita,
                       prcntOld=prcntOld,
                       percent_ba_higher=percent_ba_higher,
                       prcntWhite=prcntWhite,
                       blackpct=blackpct,
                       prcntHisp=prcntHisp,
                       people_per_mi2=people_per_mi2,
                       will_def_perc=will_def_perc,
                       will_prob_perc=will_prob_perc,
                       unsure_perc=unsure_perc,
                       not_prob_perc=not_prob_perc,
                       not_def_perc=not_def_perc
                       ) %>%
    #
    generate_predictor(time_window = seq(as.Date("2021-04-19"),unique(state_vaccines$date[state_vaccines$date %in% as.Date(lotteries$lottery_announced[lotteries$state %in% i])]),1),
                      cases = mean(cases),
                       deaths = mean(deaths)) %>%
    
    # Generate the fitted weights for the synthetic control
    generate_weights(optimization_window = seq(as.Date("2021-04-19"),unique(state_vaccines$date[state_vaccines$date %in% as.Date(lotteries$lottery_announced[lotteries$state %in% i])]),1), # time to use in the optimization task
                     optimization_method = "All", #Using all optimization functions
margin_ipop = .02,sigf_ipop = 7,bound_ipop = 6 # optimizer options
) %>%
    
    # Generate the synthetic control
    generate_control()
  
  differences <- state_vaccines_out %>% plot_differences()
  print(differences)
  
  trends <- state_vaccines_out %>% plot_trends()
  trends <- trends + labs(caption=paste("State:",i))
  print(trends)
  weights <- state_vaccines_out %>% plot_weights()
  weights <- weights + labs(caption=paste("State:",i))
  print(weights)
  balance <- state_vaccines_out %>% grab_balance_table()
  placebos <- state_vaccines_out %>% plot_placebos()
  
  ggsave(file=paste(i,"dose1_synth_differences.png",sep = "_"),differences, width = 8, height = 4.5, units = "in")
  ggsave(file=paste(i,"dose1_synth_placebos.png",sep = "_"),placebos, width = 8, height = 4.5, units = "in")
  ggsave(file=paste(i,"dose1_synth_trends.png",sep = "_"),trends, width = 8, height = 4.5, units = "in")
  ggsave(file=paste(i,"dose1_synth_weights.png",sep = "_"),weights, width = 8, height = 4.5, units = "in")


  
  administered_dose1_trends[[i]] <- trends
  administered_dose1_differences[[i]] <- differences
  administered_dose1_weights[[i]] <- weights
  administered_dose1_balance[[i]] <- balance
  administered_dose1_placebos[[i]] <- placebos
  administered_dose1_dat[[i]] <- state_vaccines_out
}

dose1_list <- list(administered_dose1_dat,
                   administered_dose1_trends,
                   administered_dose1_differences,
                   administered_dose1_weights,
                   administered_dose1_balance,
                   administered_dose1_placebos)


### Running complete series synthetic control ###

#Creating empty lists to hold results
series_complete_trends <- list()
series_complete_differences <- list()
series_complete_weights <- list()
series_complete_balance <- list()
series_complete_placebos <- list()
series_complete_dat <- list()

#setwd() #Set a working directory for images from the analyses to be saved
for(i in unique(lotteries$state[lotteries$lottery_incentive %in% 1])){
  
  # Take out all untreated states 
  
  x <- subset(state_vaccines,state_vaccines$state %in% i)
  x <- rbind(x,untreated_states)
  x <- x[order(x$date,x$state),]
  
  #take out observations before April 19th (widely available vaccine date)
  x <- x[x$date >= as.Date("2021-04-19"),]
  
  state_vaccines_out <-
    x %>%
    # initial the synthetic control object
    synthetic_control(outcome = Series_Complete_18PlusPop_Pct, # outcome
                      unit = state, # unit index in the panel data
                      time = date, # time index in the panel data
                      i_unit = i, # unit where the intervention occurred
                      i_time = unique(state_vaccines$date[state_vaccines$date %in% as.Date(lotteries$lottery_announced[lotteries$state %in% i])]), # time period when the intervention occurred
                      generate_placebos=T # generate placebo synthetic controls (for inference)
    ) %>%
    
    # Generate the aggregate predictors used to fit the weights
    generate_predictor(time_window = as.Date("2021-04-19"),
                       pres_dem_two_party_vote_percent = pres_dem_two_party_vote_percent,
                       income_per_capita=income_per_capita,
                       prcntOld=prcntOld,
                       percent_ba_higher=percent_ba_higher,
                       prcntWhite=prcntWhite,
                       blackpct=blackpct,
                       prcntHisp=prcntHisp,
                       people_per_mi2=people_per_mi2,
                       will_def_perc=will_def_perc,
                       will_prob_perc=will_prob_perc,
                       unsure_perc=unsure_perc,
                       not_prob_perc=not_prob_perc,
                       not_def_perc=not_def_perc
    ) %>%
    
    #
    generate_predictor(time_window = seq(as.Date("2021-04-19"),unique(state_vaccines$date[state_vaccines$date %in% as.Date(lotteries$lottery_announced[lotteries$state %in% i])]),1),
                       cases = mean(cases),
                       deaths = mean(deaths)) %>%
    
    # Generate the fitted weights for the synthetic control
    generate_weights(optimization_window = seq(as.Date("2021-04-19"),unique(state_vaccines$date[state_vaccines$date %in% as.Date(lotteries$lottery_announced[lotteries$state %in% i])]),1), # time to use in the optimization task
                     optimization_method = "All", #Using all optimization functions
                     margin_ipop = .02,sigf_ipop = 7,bound_ipop = 6 # optimizer options
    ) %>%
    
    # Generate the synthetic control
    generate_control()
  
  differences <- state_vaccines_out %>% plot_differences()
  print(differences)
  
  trends <- state_vaccines_out %>% plot_trends()
  trends <- trends + labs(caption=paste("State:",i))
  print(trends)
  weights <- state_vaccines_out %>% plot_weights()
  weights <- weights + labs(caption=paste("State:",i))
  print(weights)
  balance <- state_vaccines_out %>% grab_balance_table()
  placebos <- state_vaccines_out %>% plot_placebos()
  
  ggsave(file=paste(i,"complete_synth_differences.png",sep = "_"),differences, width = 8, height = 4.5, units = "in")
  ggsave(file=paste(i,"complete_synth_placebos.png",sep = "_"),placebos, width = 8, height = 4.5, units = "in")
  ggsave(file=paste(i,"complete_synth_trends.png",sep = "_"),trends, width = 8, height = 4.5, units = "in")
  ggsave(file=paste(i,"complete_synth_weights.png",sep = "_"),weights, width = 8, height = 4.5, units = "in")
  
  
  series_complete_trends[[i]] <- trends
  series_complete_differences[[i]] <- differences
  series_complete_weights[[i]] <- weights
  series_complete_balance[[i]] <- balance
  series_complete_placebos[[i]] <- placebos
  series_complete_dat[[i]] <- state_vaccines_out
}


complete_list <- list(series_complete_dat,
                    series_complete_trends,
                   series_complete_differences,
                   series_complete_weights,
                   series_complete_balance,
                   series_complete_placebos)


